

In [1]:

    
%matplotlib notebook



In [2]:

    
import numpy as np
from sympy import *
import matplotlib.pyplot as plt



In [3]:

    
init_printing()
u = Function('u')
t, eps= symbols('t epsilon')
omega = symbols('omega', positive=True)

**Note:** This notebook requires SymPy 1.5 to work.

Consider the following system

$$\ddot{u} + 4 u + \varepsilon u^2 \ddot{u} = 0 \enspace .$$

This system can be rewritten as $$(1 + \varepsilon u^2)\ddot u + 4u = 0 \enspace ,$$ or $$\ddot u + \frac{4u}{1 + \varepsilon u^2} = 0 \enspace .$$ As a first order system it reads $$\begin{pmatrix} \dot u \\ \dot v \end{pmatrix} = \begin{pmatrix} v \\ -\frac{4u}{1 + \varepsilon u^2} \end{pmatrix} \enspace ,$$ with Jacobian matrix $$J(u,v) = \begin{bmatrix} 0 & 1 \\ -\frac{4(1 - \varepsilon u^2)}{(1+ \varepsilon u^2)^2} & 0 \end{bmatrix} \enspace .$$ This system has a fixed point in $(0,0)$, with eigenvalues $$\lambda_1 = -2i,\quad \lambda_2 = 2i \enspace ,$$ and we can conclude that this fixed point is a center.



In [4]:

    
eq = (1 + eps*u(t)**2)*diff(u(t), t, 2) + omega**2*u(t)
eq









    Out[4]:





$\displaystyle \omega^{2} u{\left(t \right)} + \left(\epsilon u^{2}{\left(t \right)} + 1\right) \frac{d^{2}}{d t^{2}} u{\left(t \right)}$



In [5]:

    
ode_order(eq, u)









    Out[5]:





$\displaystyle 2$

Straightforward expansion

Let's take $u = u_0 +\varepsilon u_1 + \cdots$. Replacing this in the equation we obtain



In [6]:

    
u0 = Function('u0')
u1 = Function('u1')
subs = [(u(t), u0(t) + eps*u1(t))]



In [7]:

    
aux = eq.subs(subs)
aux.doit().expand()









    Out[7]:





$\displaystyle \epsilon^{4} \operatorname{u_{1}}^{2}{\left(t \right)} \frac{d^{2}}{d t^{2}} \operatorname{u_{1}}{\left(t \right)} + 2 \epsilon^{3} \operatorname{u_{0}}{\left(t \right)} \operatorname{u_{1}}{\left(t \right)} \frac{d^{2}}{d t^{2}} \operatorname{u_{1}}{\left(t \right)} + \epsilon^{3} \operatorname{u_{1}}^{2}{\left(t \right)} \frac{d^{2}}{d t^{2}} \operatorname{u_{0}}{\left(t \right)} + \epsilon^{2} \operatorname{u_{0}}^{2}{\left(t \right)} \frac{d^{2}}{d t^{2}} \operatorname{u_{1}}{\left(t \right)} + 2 \epsilon^{2} \operatorname{u_{0}}{\left(t \right)} \operatorname{u_{1}}{\left(t \right)} \frac{d^{2}}{d t^{2}} \operatorname{u_{0}}{\left(t \right)} + \epsilon \omega^{2} \operatorname{u_{1}}{\left(t \right)} + \epsilon \operatorname{u_{0}}^{2}{\left(t \right)} \frac{d^{2}}{d t^{2}} \operatorname{u_{0}}{\left(t \right)} + \epsilon \frac{d^{2}}{d t^{2}} \operatorname{u_{1}}{\left(t \right)} + \omega^{2} \operatorname{u_{0}}{\left(t \right)} + \frac{d^{2}}{d t^{2}} \operatorname{u_{0}}{\left(t \right)}$



In [8]:

    
poly = Poly(aux.doit(), eps)



In [9]:

    
coefs = poly.coeffs()
coefs









    Out[9]:





$\displaystyle \left[ \operatorname{u_{1}}^{2}{\left(t \right)} \frac{d^{2}}{d t^{2}} \operatorname{u_{1}}{\left(t \right)}, \  2 \operatorname{u_{0}}{\left(t \right)} \operatorname{u_{1}}{\left(t \right)} \frac{d^{2}}{d t^{2}} \operatorname{u_{1}}{\left(t \right)} + \operatorname{u_{1}}^{2}{\left(t \right)} \frac{d^{2}}{d t^{2}} \operatorname{u_{0}}{\left(t \right)}, \  \operatorname{u_{0}}^{2}{\left(t \right)} \frac{d^{2}}{d t^{2}} \operatorname{u_{1}}{\left(t \right)} + 2 \operatorname{u_{0}}{\left(t \right)} \operatorname{u_{1}}{\left(t \right)} \frac{d^{2}}{d t^{2}} \operatorname{u_{0}}{\left(t \right)}, \  \omega^{2} \operatorname{u_{1}}{\left(t \right)} + \operatorname{u_{0}}^{2}{\left(t \right)} \frac{d^{2}}{d t^{2}} \operatorname{u_{0}}{\left(t \right)} + \frac{d^{2}}{d t^{2}} \operatorname{u_{1}}{\left(t \right)}, \  \omega^{2} \operatorname{u_{0}}{\left(t \right)} + \frac{d^{2}}{d t^{2}} \operatorname{u_{0}}{\left(t \right)}\right]$



In [10]:

    
sol0 = dsolve(coefs[-1], u0(t)).rhs
sol0









    Out[10]:





$\displaystyle C_{1} \sin{\left(\omega t \right)} + C_{2} \cos{\left(\omega t \right)}$



In [11]:

    
aux = (u0(t)**2*u0(t).diff(t, 2)).subs(u0(t), sol0).doit()
aux









    Out[11]:





$\displaystyle - \omega^{2} \left(C_{1} \sin{\left(\omega t \right)} + C_{2} \cos{\left(\omega t \right)}\right)^{3}$



In [12]:

    
dsolve(u1(t).diff(t, 2) + omega**2*u1(t) + aux)









    Out[12]:





$\displaystyle \operatorname{u_{1}}{\left(t \right)} = \frac{C_{1} \left(- C_{1}^{2} + 3 C_{2}^{2}\right) \sin^{3}{\left(\omega t \right)}}{8} + \frac{3 C_{2} \left(- C_{1}^{2} + 3 C_{2}^{2}\right) \cos^{3}{\left(\omega t \right)}}{8} + \frac{C_{2} \left(3 C_{1}^{2} - 5 C_{2}^{2}\right) \cos^{5}{\left(\omega t \right)}}{8} + \left(\frac{3 C_{1}^{3}}{8} + \frac{3 C_{1}^{2} C_{2} \omega t}{8} + \frac{3 C_{1} C_{2}^{2}}{8} + \frac{3 C_{2}^{3} \omega t}{8} + C_{3}\right) \sin{\left(\omega t \right)} + \left(- \frac{3 C_{1}^{3} \omega t}{8} - \frac{3 C_{1} C_{2}^{2} \omega t}{8} - \frac{C_{2}^{3}}{2} + \frac{C_{2} \left(- 3 C_{1}^{2} + 5 C_{2}^{2}\right) \sin^{4}{\left(\omega t \right)}}{8} + C_{4}\right) \cos{\left(\omega t \right)}$



In [13]:

    
eq_aux = expand(coefs[-2].subs(u0(t), sol0))
eq_aux.doit()









    Out[13]:





$\displaystyle - C_{1}^{2} \omega^{2} \left(C_{1} \sin{\left(\omega t \right)} + C_{2} \cos{\left(\omega t \right)}\right) \sin^{2}{\left(\omega t \right)} - 2 C_{1} C_{2} \omega^{2} \left(C_{1} \sin{\left(\omega t \right)} + C_{2} \cos{\left(\omega t \right)}\right) \sin{\left(\omega t \right)} \cos{\left(\omega t \right)} - C_{2}^{2} \omega^{2} \left(C_{1} \sin{\left(\omega t \right)} + C_{2} \cos{\left(\omega t \right)}\right) \cos^{2}{\left(\omega t \right)} + \omega^{2} \operatorname{u_{1}}{\left(t \right)} + \frac{d^{2}}{d t^{2}} \operatorname{u_{1}}{\left(t \right)}$



In [14]:

    
sol1 = dsolve(eq_aux, u1(t)).rhs
sol1









    Out[14]:





$\displaystyle \frac{C_{1} \left(- C_{1}^{2} + 3 C_{2}^{2}\right) \sin^{3}{\left(\omega t \right)}}{8} + \frac{3 C_{2} \left(- C_{1}^{2} + 3 C_{2}^{2}\right) \cos^{3}{\left(\omega t \right)}}{8} + \frac{C_{2} \left(3 C_{1}^{2} - 5 C_{2}^{2}\right) \cos^{5}{\left(\omega t \right)}}{8} + \left(\frac{3 C_{1}^{3}}{8} + \frac{3 C_{1}^{2} C_{2} \omega t}{8} + \frac{3 C_{1} C_{2}^{2}}{8} + \frac{3 C_{2}^{3} \omega t}{8} + C_{3}\right) \sin{\left(\omega t \right)} + \left(- \frac{3 C_{1}^{3} \omega t}{8} - \frac{3 C_{1} C_{2}^{2} \omega t}{8} - \frac{C_{2}^{3}}{2} + \frac{C_{2} \left(- 3 C_{1}^{2} + 5 C_{2}^{2}\right) \sin^{4}{\left(\omega t \right)}}{8} + C_{4}\right) \cos{\left(\omega t \right)}$



In [15]:

    
C1, C2, C3, C4 = symbols("C1 C2 C3 C4")



In [16]:

    
print(sol1.subs({C1: 0, C2: 1, C3: 0, C4: -S(1)/4}))









    



3*omega*t*sin(omega*t)/8 + (5*sin(omega*t)**4/8 - 3/4)*cos(omega*t) - 5*cos(omega*t)**5/8 + 9*cos(omega*t)**3/8



In [ ]:



In [17]:

    
u_app = sol0 + eps*sol1



In [18]:

    
aux_eqs = [
        sol0.subs(t, 0)-1,
        sol1.subs(t, 0),
        diff(sol0, t).subs(t, 0),
        diff(sol1, t).subs(t, 0)]
aux_eqs









    Out[18]:





$\displaystyle \left[ C_{2} - 1, \  - \frac{C_{2}^{3}}{2} + \frac{3 C_{2} \left(- C_{1}^{2} + 3 C_{2}^{2}\right)}{8} + \frac{C_{2} \left(3 C_{1}^{2} - 5 C_{2}^{2}\right)}{8} + C_{4}, \  C_{1} \omega, \  - \frac{3 C_{1}^{3} \omega}{8} - \frac{3 C_{1} C_{2}^{2} \omega}{8} + \omega \left(\frac{3 C_{1}^{3}}{8} + \frac{3 C_{1} C_{2}^{2}}{8} + C_{3}\right)\right]$



In [19]:

    
coef = u_app.free_symbols - eq.free_symbols
coef









    Out[19]:





$\displaystyle \left\{C_{1}, C_{2}, C_{3}, C_{4}\right\}$



In [20]:

    
subs_sol = solve(aux_eqs, coef)
subs_sol









    Out[20]:





$\displaystyle \left\{ C_{1} : 0, \  C_{2} : 1, \  C_{3} : 0, \  C_{4} : 0\right\}$



In [21]:

    
u_app2 = u_app.subs(subs_sol)
trigsimp(u_app2)









    Out[21]:





$\displaystyle \epsilon \left(\frac{3 \omega t \sin{\left(\omega t \right)}}{8} + \frac{\left(5 \sin^{4}{\left(\omega t \right)} - 4\right) \cos{\left(\omega t \right)}}{8} - \frac{5 \cos^{5}{\left(\omega t \right)}}{8} + \frac{9 \cos^{3}{\left(\omega t \right)}}{8}\right) + \cos{\left(\omega t \right)}$



In [22]:

    
print(u_app2)









    



epsilon*(3*omega*t*sin(omega*t)/8 + (5*sin(omega*t)**4/8 - 1/2)*cos(omega*t) - 5*cos(omega*t)**5/8 + 9*cos(omega*t)**3/8) + cos(omega*t)



In [23]:

    
final_sol = trigsimp(u_app2).subs(omega, 2).expand()
final_sol









    Out[23]:





$\displaystyle \frac{3 \epsilon t \sin{\left(2 t \right)}}{4} + \frac{5 \epsilon \sin^{4}{\left(2 t \right)} \cos{\left(2 t \right)}}{8} - \frac{5 \epsilon \cos^{5}{\left(2 t \right)}}{8} + \frac{9 \epsilon \cos^{3}{\left(2 t \right)}}{8} - \frac{\epsilon \cos{\left(2 t \right)}}{2} + \cos{\left(2 t \right)}$



In [24]:

    
trigsimp(final_sol).expand()









    Out[24]:





$\displaystyle \frac{3 \epsilon t \sin{\left(2 t \right)}}{4} - \frac{\epsilon \cos^{3}{\left(2 t \right)}}{8} + \frac{\epsilon \cos{\left(2 t \right)}}{8} + \cos{\left(2 t \right)}$



In [25]:

    
sol = (1 - 5*eps/32)*cos(2*t) + eps/32*(6*cos(2*t) - cos(6*t) + 24*t*sin(2*t))
sol









    Out[25]:





$\displaystyle \frac{\epsilon \left(24 t \sin{\left(2 t \right)} + 6 \cos{\left(2 t \right)} - \cos{\left(6 t \right)}\right)}{32} + \left(1 - \frac{5 \epsilon}{32}\right) \cos{\left(2 t \right)}$



In [26]:

    
plot((sol - final_sol).subs(eps, 1))









    














    











    Out[26]:





<sympy.plotting.plot.Plot at 0x7f85a59ae4c0>



In [27]:

    
from scipy.integrate import odeint
def fun(x, t=0, eps=0.1):
    x0, x1 = x
    return [x1, -4*x0/(1 + eps*x0**2)]

t_vec = np.linspace(0, 100,  1000)
x = odeint(fun, [1, 0], t_vec, args=(0.1,))



In [28]:

    
lam_sol = lambdify((t, eps), final_sol, "numpy")
uu = lam_sol(t_vec, 0.1)



In [30]:

    
plt.figure()
plt.plot(t_vec, x[:,0])
plt.plot(t_vec, uu, '--r')
plt.show()



In [ ]: